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Abstract 

We propose a new formulation for integrating over smooth curves and surfaces 
that are described by their closest point mappings. Our method is designed for 
curves and surfaces that are not defined by any explicit parameterization and is 
intended to be used in combination with level set techniques. However, contrary 
to the common practice with level set methods, the volume integrals derived 
from our formulation coincide exactly with the surface or line integrals that one 
wishes to compute. We study various aspects of this formulation and provide a 
geometric interpretation of this formulation in terms of the singular values of the 
Jacobian matrix of the closest point mapping. Additionally, we extend the 
formulation - initially derived to integrate over manifolds of codimension one - to 
include integration along curves in three dimensions. Some numerical examples 
using very simple discretizations are presented to demonstrate the efFicacy of the 
formulation. 

Keywords: boundary integrals; closest point mapping; level set methods 


1 Introduction 

This paper provides simple formulations for integrating over manifolds of codimen¬ 
sions one, or two in when the manifolds are described by functions that map 
points in R" (n = 2,3) to their closest points on curves or surfaces using the Eu¬ 
clidean distance. The idea for the present work originated in [1] where the authors 
proposed a formulation for computing integrals of the form 

f v{x{s))ds, (1) 

JOQ 

in the level set framework, namely when the domain is represented implicitly by 
the signed distance function to its boundary dfl. Typically in a level set method 
[2, 3, 4], to evaluate an integral of the form of (1) where dfl is the zero level set 
of a continuous function ip, it is necessary to extend the function v defined on the 
boundary to a neighborhood in R”. The extension of v, denoted v, is typically a 
constant extension of v. The integral is then approximated by an integral involving 
a regularized Dirac-i5 function concentrated on 90, namely 




i;(x)(5e((p(x))|V(p(x)|dx. 


Various numerical approximations of this delta function have been proposed, see 
e.g. [5, 6, 7, 8, 9], 
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In [1], with the choice of (/? = dga being a signed distance function to did, the 
integral ( 1 ) is expressed as an average of integrals over nearby level sets of dgn, where 
these nearby level sets continuously sweep a thin tubular neighborhood around the 
boundary Oil of radius e. Consequently, (1) is equivalent to the volume integral 
shown on the right hand side below: 

[ v{x{s))ds = [ i;(x*)J(x;dan)^e(dan(x))dx, ( 2 ) 

Jan Jm" 

where 6^ is an averaging kernel, x* is the closest point on dil to x and J(x;daa) 
accounts for the change in curvature between the nearby level sets and the zero 
level set. 

Now suppose that dil is a smooth hypersurface in and assume that x is suffi¬ 
ciently close to n so that the closest point mapping 


X* = fbn(x) = argmin^gaolx - y| 


is continuously differentiable. Then the restriction of PgQ to dilr) is a diffeormor- 
phism between dilri and dil, where dilri ■= {x : dgQ{x.) = 77 }. As a result, it is 
possible to write integrals over dil using points on dilri as: 


f v{x)dS = f 
J dCl J dQri 


u(x*) J(x; ri)dS, 


where J (x, 77 ) comes from the change of variable defined by PgQ restricted on dilri ■ 
Averaging the above integrals respectively with a kernel, Se, compactly supported 
in [—e, e], we obtain 


f v{x.)dS = f Se(ji) f v{x*)J(x;r])dS dr], 
J dO, J — e J dQ,Ti 


Formula (2) then follows from the coarea formula [10] applied to the integral on the 
right hand side. 

In the following section, we show that in three dimensions the Jacobian J in (2) 
is the product of the first two singular values, cti and a 2 , of the Jacobian matrix of 
the closest point mapping namely. 


f vix{s))ds = f 
JdQ Jr- 


2 

v{Pgn{x))S^{dgn{x)) ]^crj(x)dx. 

i=i 


(3) 


To motivate the new approach using singular values, we consider Cartesian coor¬ 
dinate systems with the origin placed on points sufficiently close to the surface, 
and the 2 : direction normal to the surface. Thus the partial derivatives of the clos¬ 
est point mapping in the z direction will yield zero and the partial derivatives in 
the other two directions naturally correspond to differentiation in the tangential 
directions. Thus we see that one of the singular values should be 0 while the other 
two are related to the surface area element. We also derive a similar formula for 
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integration along curves in three dimensions (codimension 2). The advantages of 
this new formula include the ease for constructing higher order approximations of 
J via e.g. simple differencing, even in neighborhoods of surface boundaries where 
curvatures become unbounded. 

This paper is motivated by the recent success in the closest point methods and 
the Dynamic Surface Extension method [11], for evolving interfaces and solving 
partial differential equations on surfaces [12, 13, 14, 15], by the need to process 
data sets that contain unstructured points sampled from some underlying surfaces, 
and targets applications where manifolds are not defined by patches of explicit 
parameterizations and may evolve drastically due to some coupled processes; see 
e.g. free boundary problems [16]. Our work provides a convenient way to formu¬ 
late boundary integral methods in such applications without conversion to local 
parameterizations. If the manifolds are defined by explicit parameterizations, it is 
natural and typically more accurate to use conventional methods such as Nystrom 
methods using quadratures on the parameter space or Boundary Element Methods 
with weak formulations, see e.g. [17]. Additionally, for applications involving fluid- 
structure interactions, we mention the immersed boundary method which involves 
accurate discretizations of Dirac delta measures [18, 19]. 

Closest point mappings are easily computed in the context of level set methods 
[3] since there exist fast algorithms for constructing distance functions from level 
set functions [20, 21, 22, 23, 24]. More precisely, 

Pan{^) = X - (ian(x)Vdan(x). 

Our previous work [1] as well as this current paper provide a simple framework for 
constructing numerical schemes for boundary integral methods when the interface is 
described implicitly by a level set function, and is intended for use in such context. 

Finally, closest point mappings can also be computed easily from dense and unor¬ 
ganized point sets that are acquired directly from an imaging device (e.g. LIDAR). 
This paper lays the foundation of a numerical scheme for computing integrals over 
surfaces sampled by unstructured point clouds. 

2 Integration using the closest point mapping 

In this section, we relate the Jacobian J in (2) to the singular values of the Jacobian 
matrix of the closest point mapping from or to T, where T denotes the curves 
or surfaces on which integrals are defined. We assume that in three dimensions, if 
r is not closed, it has smooth boundaries. For clarity of the exposition in the rest 
of the paper, we will now denote the distance function simply by d. 

2.1 Codimension 1 

We consider a compact curve or surface T that can either be closed or not. If T 
is closed, then it is the boundary of a domain D so that T can be denoted dO,. If T is 
not closed, we assume that it has smooth boundaries. We define d : R" i—>• RU{0} to 
be the distance function to T and Pr to be the closest point mapping Pr : R" i—J- T 
(for n = 2,3) defined as 


|Pr(x) -x| =min|y-x|. 
yer 


( 4 ) 
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We let do be the distance function to F if it is open and dg be the signed distance 
function to F = 917 if it is closed. The signed distance function is defined as 


ds(x) 


infygoc |x - y| if x G 17, 
- infygn |x - y| if x G 17^. 


Then we define d as follows: 


d(x) 


do (x) if F is open, 
ds{x) if F is closed. 


(5) 


The following lemma provides a concise expression of the Gaussian curvature in 
terms of the distance function. This is probably a known result but we include its 
proof to preserve the completeness of the paper. 


Lemma 1 Let d be the distance function to F defined in (5). For \r]\ sufficiently 
close to 0, the Gaussian curvature at a point on the rj level set F,, := : d{f) = rj} 

can be expressed as 


Gij dxxdyy dxxdzz “f dyydz^z ^xy ^xz ^yz' 


( 6 ) 


Proof Starting with the definition of the Gaussian curvature G for a surface (see 
[25]), we can obtain an expression for the Gaussian curvature of its 77 -level set in 
terms of d as 


G = {'Vd,adj{Hess{d))Vd) 

~ d-xi^yy^zz — dyfi) + dy{dxxdzz — dxfij + dfidxxdyy — dxy) 

T ‘2fixdy(^dxzdyz dxydzz) F dydzi^dxydxz dyzdxx^ 

T dxdz^dxydyz dxzdyy^f (7) 


We show that this expression is the same as ( 6 ) by rearranging the terms above 
and using the fact that close to F the distance function satisfies |Vd| = 1. First we 
rearrange the terms in G: 




dydx 


dzz + d^dxxdyy dxdy^ dydxz d^dxy 
F^\dxdy{^dxzdyz dxydzz^ F dydzi^dxydxz dyzdxx^ F dxdzi^dxydyz dxzdyy)f 


and rewrite each of the first six terms in terms of jVdp, e.g. 


^x^yy^zz 


= \Vd\^d 

yy^ZZ dydyydzz d^dyyd^Z dyyd^^ dydyyd^Z d^dyyd^Z ’ 


= 1 
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Thus we have 


z^xy 


^x^yy^zz T dydxx^zz T d^dxxdyy d^d^,^ d^,d^^ d^d, 

dxxdyy T dxxdz 


-*^xx'~^yy 


^x'^yz 

72 ^2 

-^xz 


'^y^xz 

72 . 


dyydzz d^y d^^ dy^ dydyydzz d^dyydzz 


— Gy 

d^dxxdzz d^dxxdzz dydxxdyy dj.dxxdyy 


( 8 ) 


+dydyz + d^dy^ + d^d^^ 


+ d^zdi, + did: 


‘X ^xy 


d^d^ 

^y^xy 


Using ( 8 ) and rearranging the rest of the terms in (7) we obtain G = Gy. □ 

Proposition 2 Consider a G^ compact surface F C M" (n = 2,3) of codimen¬ 
sion 1 and let d be defined as in (5). Define the closest point projection map Pr as 
in (4) for x G R". For \p\ sufficiently close to zero, let F^ be the rj level set of d 


Ty := {x : d(x) = 77 } . 


Define the Jacobian Jy as 


Jy . - 


1 + pKy if n = 2, 

l-\-2rjFIy + rj'^Gy if n = 3, 


(9) 


where Ky is the signed curvature ofTy in 2D, and Hy and Gy are its Mean curvature 
and Gaussian curvature respectively in 3D. 

Then if Pf is the Jacobian matrix of Pr we have 


Jn 


ai, n = 2, 
(71(72, n = 3, 


( 10 ) 


where ai, are the first two singular values of the Jacobian matrix Pf. 


Proof The distance function d satisfies the property d{x.) = 0 for x G F. Also, since 
F is C^, its distance function d belongs to see e.g. [26, 27]. It follows 

that the order of the mixed partial derivatives does not matter. In addition, the 
normals to a smooth interface do not focus right away so that the distance function 
is smooth in a tubular neighborhood T around F, and is linear with slope one along 
the normals. Therefore we have 


jVdj = 1 for all x G T. (II) 

The third important fact is that the Laplacian of d at a point x gives (up to a 
constant related to the dimension) the mean curvature of the isosurface of d passing 
through X, namely 


Ad(x) = (I — n)iF(x), 


( 12 ) 
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where H{x.) is the Mean curvature of the level set {y : d{y) = d{x.)}. Differentiating 
(11) with respect to each variable gives the following equations in three dimensions: 

dxdxx dydxy “t“ dzdxz — 0, (13) 

dxdyx dydyy “h dzdyz — 0, 

dxdzx dydzy ~\~ dzdzz — 0 - ( 1 ^) 

In particular the two dimensional case can be derived by assuming that the distance 
function is constant in z. 


Two dimensions. In that case the Jacobian matrix Pp of the closest point pro¬ 
jection map is 

p/ _ f 1 dj. ddxx (^dydx H“ ddyx^ j 

y i^dxdy -|- ddxy^ 1 dy ddyy J 

Since Schwartz’ Theorem holds, we have dxy = dyx making Pp a real symmetric 
matrix. It is therefore diagonalizable with eigenvalues 0 and 1 — dAd. Indeed, we 
have 


PAVd = 


^ dxi, 1 dj. dy ) d( dxdxx T dydyx ) 

-0 by ( 11 ) in 2D -o by (is) in 2D 

dy ( 1 d^ dy ) d( dydyy “h dxdxy ) 

y ^0 by ( 11 ) in 2D -o by (i4) in 2D ) 


= 0 , 


and for v = 


dll 


PpV = 


dy T dyd^ dyddxx dj.dy ddxdxy 

dydx T dyddxx dydx dxddyy 


dy 


dydxx dxdxy 

dydxy dxdyy 


= V 


^ Ad{ dy) (dydyy “h dxdxy^ 
^0 by (14) in 2D 

Ad(dx) “h dxdxx T dydxy 


V 

= (1 — dAd)-v. 


=0 by (13) in 2D J 


Since ||v|| = 1, v is an eigenvector corresponding to the eigenvalue A = 1 — dAd. 
Thus, for X such that d(x) = rj we have that the eigenvalue A of Pp satisfies 


A = 1 — rjAd = 1 + rjKrj 
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by (12). Since 1 + rjKrj > 0, it follows that A coincides with the singular value of Pp 
and hence 


tTl = 1 + TJKn- 


Three dimensions. Since for \r]\ sufficiently close to 0 the distance function is C*^, 
the Jacobian matrix 



Pp (^dxdy ddxy^ 1 dy ddyy {dzdy ddzy^ 

y (.dxdz ddxz^ “t“ ddyz^ 1 d^ ddzz 


is a real symmetric matrix which is diagonalizable with one zero eigenvalue and two 
other eigenvalues Ai and A 2 . Indeed using (13),(14),(15) and (11) we can show that 


PpVd = 0. 


Now consider x such that d(x) = rj. Then, the characteristic polynomial x(A) of Pp 
is 


X(A) = -A (A^ - (2 - ryAd)A - Q ), 


where Q = —Gyrf’+rjlS.d— 1 with Gy defined in ( 6 ). Since the other two eigenvalues 


of Pp are the solutions of the quadratic equation A^ — (2 — riAd)X — Q = 0, it 
follows that 

A1A2 = -Q = 1 - rjAd + rj'^Gy = 1 + 2 r]Hy + if Gy. 

Since 1 + 2riHy + rf'Gy > 0, it follows that 
tTicr2 = 1 + 2r]Hy + rfGy, 

where ui and 02 are singular values of Pp. □ 

This leads to the following proposition: 

Theorem 3 Consider T a curve in 2D or surface in 3D with G^ boundaries if it 
is not closed, and define d : M" 1 —R'*" U {0} (n = 2,3) to be the distance function to 
r with Pr : R” 1 —>■ T the closest point mapping to T. Then for emaxa;gr |k( 2 :)| < 1 
for any nfx) principal curvatures of T at x, we have 



( 16 ) 


where 6^ is an averaging kernel and S(x)is defined as 



cri(x), n = 2 , 

(Ti(x)(T 2 (x), n = 3. 
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where o'j(x) , j = 1, 2, is the j-th singular value of the Jacobian matrix Pf evaluated 
at X. 

Proof If r is closed we combine Equation (2) with the result J(x) = E(x) from 
Equation (10) of Proposition 2. 

If r is open there is a little more to show since Equation (2) was only derived for 
closed manifolds. Before we state the result, it is necessary to understand how E^ 
defined in (9) (an 77 —level set of d) looks like for an open curve in two dimensions 
and for a surface with boundaries in three dimensions. 

In two dimensions, E^ consists of a flat tubular part on either side of the curve and 
two semi circles at the two ends of the curve. See Eigure 1 . 

In three dimensions E is in general made up of three distinct parts: the interior 
part, the edges of the boundary and the corners. If we assume that E has N edges 
then we can write E = E° U U where E° is the interior of E, Ei 

is the f-th edge of the boundary of E and Ci is its i -th corner. In that setting we 
can write E^ = /,, U U (u(^;^S'f), where In is the inside portion of E^, Tf 

is the cylindrical part of E^ representing the set of points located at a distance ry 
from the i -th edge Ei, and finally is the spherical part of E^ representing the 
set of points located at a distance ij from the i -th corner Ci. See Eigure 2. 

In both cases we need to integrate over E^ and then subtract the two semi circles 
at the two end points of the curve (in two dimensions) or subtract the portions 
of sphere at the corners of the surface and the portions of cylinders at the edges 
of the surface (in three dimensions). However, it turns out that the subtraction is 
unnecessary since E(x) = 0 on each of the subtracted pieces as shown below. 

Two dimensions. On the semi-circle around the end point of a curve, the closest 
point mapping is constant since all points on the semi-circle E^ map to the end 
point. As a result, the singular values of the Jacobian matrix of the closest point 
mapping are all zeros and thus S(x) = 0 on the semi-circles around the end points 
of a curve. 

Three dimensions. As in two dimensions, on the portions of sphere around a cor¬ 
ner point of a surface, the closest point mapping is constant and thus S(x) = 0. 
On the portion of cylinders, the closest point mapping is constant along the ra¬ 
dial dimension (one of the principal directions or singular vector) resulting of the 
singular value along that direction to be zero. Since E(x) is the product of the 
singular values, it follows that E(x) = 0 on the portion of cylinders as well. Con¬ 
sequently, Equation (16) holds for any curve or surface with boundaries of 
codimension 1. □ 

2.2 Codimension 2 

We consider a curve in denoted by E and let 7 ( 3 ) be a parameterization by 
arclength of E. We denote by d : 1 —>■ R"*" U {0} the distance function to E and let 
Pr : E be the closest point mapping to E. We consider a parameterization of 

the tubular part of the level surface for rj G [ 0 , e] defined as 

x(s, 9, 7 ) : 7 ( 5 ) -I- r]COs6'N{s) + rj sin 9'B{s), 
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where T = N and B constitute the Frenet frame for 7 as illustrated in Figure 3. 
As in the previous section, if F is closed then d is the signed distance function to F. 

If we project a point x on the tubular part of the level surface F^ defined in (9), 
we have iA(x(s, 9,r])) = 7 (s). If L is the length of the curve it follows that 


277 pL 

0 Jo 


nlTT pL 

5 (Pr(x(s, 0 , ? 7 )))|xs X xe|dsd 6 » = / / g{s))ri{1 - r]K{s) cos 0)dsde, 

Jo Jo 


(17) 


pL/ p ZTi 

V / / {1 — r]K cos 9)d9ds 

Jo Jo 

j 9 his))ds. 


= 2717] 


Note that the tubular part of the level surface F^ does not contain the two hemi¬ 
spheres of F^ which are located at the two end points of the curve F. Thus, 


r,\{CiUC 2 } 


5(Pr(x))d5x 


2777] 



(18) 


where Ci and C 2 are the two hemispheres of the level surface F^ located at the two 
end points of the curve F. Consequently, for sufficiently small e and by the coarea 
formula we obtain 


[ 9{lis))ds = ^ f \~ [ 9{Pr{x.))] Keiv)dv, 

Jr 2,77 Jq \V Jr^\{CiUC 2 } J 

= ^ [ ff(-Pr(x))^^^^X(CiuC2)‘=(x)dx, 

Z77 d 

where is a averaging kernel supported in [0,e] and X(CiUC 2 )‘'(^) char¬ 
acteristic function of the set (Ci U C 2 Y- Because of the term in the above 

equation and for better accuracy, we choose a kernel that satisfies the condition 
iF((0) = 0. In our numerical simulations we consider the kernel 

K^iv) = “ (1 “ ‘= 0 ® (^^e)) (19) 

Since the formulation above does not use the two hemispheres located at both 
end points of the curve, in order to integrate over the tubular part of F^ only, it is 
necessary to subtract the integration over each of the hemispheres Ci and C 2 ■ The 
result can be summarized in the following proposition: 

Proposition 4 Consider a single curve F in K® parameterized by 7 ( 3 ) where 
s is the arclength parameter, and let d be the distance function to F. We define 
to be a averaging kernel compactly supported in [0, e] and Pp : '—!> B to be the 

closest point mapping to F. 
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If g is a continuous function defined on T then for sufficiently small e > 0 we 
have 


( 9 i.li.s))ds = f j g{x^)riK^{r))dT], ( 20 ) 

Jr 27r Jr 3 d(x) Jq 

where x,, is a point on a sphere of radius rj. 

Note that for the computation of the length of a curve, the correction terms given 
by integrating over both Ci and C 2 is 

rtsM4v<i,=2,. 

Jo V Jo V 

This simple correction is, however, not suitable for more general cases that contain 
multiple curve segments and several integrands. We shall derive a more elegant and 
seamless way to perform such correction in the following section. 

Now if we consider a curve in three dimensions and let Pp be its closest point 
mapping, we have the following proposition: 

Theorem 5 Let cr(x) he the nonzero singular value of PJ and let g be a contin¬ 
uous function defined on T. If 7 ( 3 ) is the arclength parameterization of T and if 
emaxa;gr |'«(a^)| < 1; where k { x ) is the curvature of the curve at x, we have 

9 il{s))ds = g(Pr(x))^^^cr(x)dx, (21) 

where d is the distance function to T. 

Proof Since is compactly supported in [0, e] it is sufficient to consider points in 
the tubular neighborhood of the curve T. Thus, for x in the tubular neighborhood, 
there exists 0 < ry < e such that x G T^. 

Case 1: x is on the spherical part of T^ corresponding to the ry-distance to either 
of the two end points of the curve T. WLOG we assume that x is at a distance g 
from the first end point Ci parameterized by 7 ( 0 ). The result is the same if x is 
on the other sphere, i.e. at a distance ry from the other end point C 2 - In that case, 
Pr(x) = 7 ( 0 ) for all x on the spherical part so that the Jacobian matrix Pf = 0. 
Therefore, for x on the spherical part of T^, all singular values of the Jacobian 
matrix are zero. 

Case 2: x is on the tubular part of T^. In that case, if we use the Frenet frame 
centered at the point x = x(s, 0 ,ry) G T^ , we can write x in the new coordinate 
system (T,N, B) as 


X = 7 ( 5 ) + uN + wB, 


( 22 ) 
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where m = 0 is the coordinate of x along T, n is the coordinate along N and w is 
the coordinate along B. Since the projection -Pr(x) = 7 (s) does not depend on v 
nor w (since the plane (N,B) is normal to the curve F) it follows that 


gPr(x) ^ 9Pr(x) 
dv dw 

On the other hand, we have 


9Pr(x) d^{s) ds dj{s) ds ^ 

du du du ds du ’ 

where ^ is the variation of the arclength parameter s with respect to u when the 
point X is moving on F,; along the tangential direction T. Since u is the arclength 
parameter along the tangential direction T, it follows that we have a unit speed 
parameterization along T giving the identity 



1 . 


In addition, 


ax dj(s) N B 
^ “ ds 

= T — kuT + ruB — rruN 
= (1 — Kv) T — rwN + ruB, 


where k is the curvature of F at 7 ( 5 ) and r is the torsion of the curve F at the 
point 7 ( 3 ). Since the level surface F,; is a tube of radius 7 , its intersection with the 
normal plane (N, B) is a circle of radius rj. Hence if we use polar coordinates on 
the normal plane, we obtain v = rj cos 9 and w = rj sin 6 . It follows that 


ax ^ 

— -T: = 1 — nrj cos 0. 
os 


Consequently we have 



ds dx ^ 
du ds 


du 


(1 — Krj cos 9), 


and 


^ _ I 
du 1 — Krj cos 9 


Therefore, in the Frenet frame, the Jacobian matrix of the closest point projection 
map can be written as 


Pr = 


/ 

V 


1 

l — KT] COS 9 
0 
0 


0 

0 

0 


0 

0 

0 
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where is the nonzero eigenvalue of the Jacobian of the closest point map¬ 

ping. Based on the hypothesis on the size of e related to the geometry of the curve 
r, the term is strictly positive and therefore is also the singular value cr(x) 

of the Jacobian of the closest point mapping. 

Therefore we have 

0 

1 

l — KT) COS B 



if X is on the spherical part of , 
if X is on the tubular part of . 


(23) 


Now using (17) and (18) we obtain 


[ 5(-fr(x))cr(x)dS'x 
Jr„ 


'r,\{CiUC2} 
o27t pL 




nZ'K nLi 

/ / g(Pr(x))cr(x)|xs X xg|dsd6> 

Jo Jo 


1 — r7K(s) COS0 


dsdO 


= 27rry / g{j{s))ds 


It follows that for a averaging kernel compactly supported in [0,e], for suffi¬ 
ciently small e and by the coarea formula, we have 



^ f - f 9iPr{x))a{x)K4'q)dr] 
271 Jo V Jrr, 

[ 5(^r(x))^^^^CT(x)dx. 
ztt Jr3 a 


□ 


3 Numerical simulations 

In this section we investigate the convergence of our numerical integration using 
simple Riemann sums over uniform Cartesian grids. Unless stated otherwise, the 
singular values are computed from the matrix the elements of which are computed 
by the standard central difference approximations of the Jacobian matrix Pp. In 
other words, the Jacobian matrix Pp is computed by using finite differences to 
evaluate the partial derivatives of each component of Pp(x); more precisely, if 
Pr(x) = (pi(x),p 2 (x),p 3 (x)), and x = (a;i,X 2 ,a; 3 ) we use finite difference to 

approximate 7 -^ for 1 < j,k < 3. We do not evaluate the expressions that involve 
OXk 

the partial derivatives of the distance function. 

In our computations we use the cosine kernel 

Ke^^iv) = Xl-e,e]{v)^ (1 + cos (^)) (24) 

for integration on surfaces of codimension 1, and the kernel 4 defined in (19) for 
codimension 2. With these compactly supported kernels, formulas (16) and (21) can 







Kublik and Tsai 


Page 13 of 17 


be considered integration of functions defined on suitable hypercubes, periodically 
extended. In such settings, simple Riemann sums on Cartesian grids are equivalent 
to sums using Trapezoidal rule, and if all the terms are known analytically, the 
order of accuracy will be related in general to the smoothness of the integrands; 
exception can be found when the normals of the surfaces are rationally dependent 
on the step sizes used in the Cartesian grids. 

3.1 Integration of codimension one surfaces 

We tested our numerical integration on two different portions of circle, a torus, a 
quarter sphere and a three quarter sphere. We computed their respective lengths or 
surface areas by integrating the constant 1 over the curve or surface. Each of these 
tests were designed to exhibit the convergence rate of our formulations on cases with 
varying difficulty. In particular, the convergence rate of our formulation depends on 
the smoothness of the closest point mapping inside the tubular neighborhood of the 
curve or surface. 

The results for the portions of circle are given in Tables 1 and 2. In the first 
convergence studies (Table 1) the line where the closest point mapping has a jump 
discontinuity is parallel to the grid lines. In this case we see a second order con¬ 
vergence rate using central differencing to compute the Jacobian matrix P^. In the 
second test case however, the portion of circle is chosen so that the line where the 
closest point mapping has a jump discontinuity is not parallel to the grid lines. In 
that case the normal to the curve is rationally dependent on the step size of the 
Cartesian grid and the convergence rate reduces to first order even though we used 
central differencing to compute Pp. We note that in these two tests, we chose e (the 
half width of the tubular neighborhood around the curve) small enough so that the 
line where the closest point mapping is discontinuous is outside of it. 

In three dimensions we first tested our method on a torus (closed smooth surface). 
The results for the torus are reported in Table 3. In this case the closest point 
mapping is very smooth and we see third order convergence when using the exact 
signed distance function and a third order difference scheme to approximate Pp 
(see REqo in Table 3). We also tested our method with a computed signed distance 
function. We constructed the signed distance function using the algorithm described 
in [ 20 ] , and compared the performance of our method with a fourth order accurate 
signed distance function and a first order accurate signed distance function (see RE 4 
and REi in Table 3.) With the fourth order accurate signed distance function we 
used a third order accurate difference scheme to approximate Pp, and with the first 
order accurate signed distance function we used a second order accurate difference 
scheme to approximate Pp. 

For surfaces with boundaries we tested the method on a quarter sphere and a 
three quarter sphere. The three quarter sphere case is illustrated in Figure 4. The 
reason for choosing these two cases is because the closest point mapping has a 
different degree of smoothness for each of these surfaces. For the quarter sphere 
the closest point mapping is smooth enough, but for the three quarter sphere, the 
tubular neighborhood around the surface contains the line where the closest point 
mapping has a jump discontinuity. In that latter case, it is therefore necessary to 
use an adequate one sided discretization to compute Pp accurately. The one-sided 
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discretization that we used is reported in Section 3.3. The test for the quarter sphere 
still uses central differencing to compute Pp. The results for the portions of sphere 
are reported in Tables 4 and 5. 

3.2 Integrating along curves in three dimensions 

In codimension 2, we tested our numerical integration on a coil wrapped around 
the helix defined parametrically as 

x(t) = {r cos{t),r sm{t),bt), 

with r = 0.75 and b = 0.25. The coil is then wrapped around the helix at a distance 
of 0.2 from the helix. See Figure 5. As our test case, we computed the length of the 
coil by integrating 1 along the curve. The results are reported in Table 6. 

3.3 One-sided discretization of the Jacobian matrix 

Here for completeness, we describe the one-sided discretization used in computing 
results reported in Table 5. For simplicity we provide the explanation in The 
discretization generalizes easily to 3D. 

We will describe the one-sided discretization for a uniform Cartesian grid in 
namely for Pr(xi j) = (Uij^Vij) with x^j = (ih,jh), i,j € Z and h > 0 being 
the step size. The Jacobian matrix will be approximated by simple finite differences 
defined below: 



The discretization of U and V have to be defined together because the two functions 
are not independent of each other. With 



and the smoothness indicator 



we define 


{Ux)i,j 



(t/j, )ij-, otherwise, 


and (Vx)ij is defined according to the choice of stencil based on S^{Uij) 




(Kr )i,jj otherwise. 


The discretization of Uy and Vy is defined similarly with the choice of the stencil 
based on S^{Vij). 
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4 Summary 

In this paper, we presented a new approach for computing integrals along curves 
and surfaces that are defined either implicitly by the distance function to these 
manifolds or by the closest point mappings. We are motivated by the abundance 
of discrete point sets sampled from surfaces using devices such as LIDAR, the 
need to compute functionals defined over the underlying surfaces, as well as many 
applications involving the level set method or the use of closest point methods. 

Contrary to most other existing approximations using either smeared out Dirac 
delta functions or locally obtained parameterized patches, we derive a volume inte¬ 
gral in the embedding Euclidean space which is equivalent to the desired surface or 
line integrals. This allows for easy construction of higher order numerical approx¬ 
imations of these integrals. The key components of this new approach include the 
use of singular values of the Jacobian matrix of the closest point mapping, which 
can be computed easily to high order even by simple finite differences. 
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Figures 


Figure 1 Level set of a 2D open curve. An example of an open curve F (black curve) and its 
7 /-level set Fr] (red curve). consists of a tubular part and two semi circles at the two ends. 


Figure 2 Level set of a 3D surface with boundaries. An example of a surface with boundaries 
viewed from different angles and its corresponding 77 -level set Frj viewed from the same angles. 
The figure at the bottom right corner shows the surface and Frj. 


Figure 3 Level set of an open curve in 3D. Three dimensional curve with its 77 -level surface Frj in 
green and the Frenet frame at a point on Frj. 


Figure 4 Three quarter sphere. The three quarter sphere and its corresponding 77 -level set F^. 


Figure 5 Coil and one of its level sets. The coil and one of the level sets of the distance function 
to the coil used in the reported numerical simulations. 


Tables 

Table 1 Errors for a portion of circle. Relative errors in the numerical approximation of the length of 
a planar curve, which is a portion of circle of radius R = 0.75 centered at 0. The width for the 
tubular neighborhood of the curve is e = 0.2. In this computation, the closest point mapping has a 
jump discontinuity along a straight-line which is arranged to be parallel to the grid lines. 


n 

Relative Error 

Order 

64 

2.7994 X 10-'‘ 

- 

128 

7.0665 X lO-'" 

1.99 

256 

1.7187 X 10““ 

2.04 

512 

4.2719 X 10-'= 

2.01 

1024 

1.0636 X 10“'= 

2.01 

2048 

2.6567 X 10“' 

2.00 

4096 

6.6045 X 10““ 

2.01 

8192 

1.6513 X 10““ 

2.00 
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Table 2 Errors for a tilted portion of circle. Relative errors in the numerical approximation of the 
length of a planar curve, which is a portion of circle of radius R = 0.75 centered at 0. The width for 
the tubular neighborhood of the curve is e = 0.2. In this computation, the jump discontinuity of the 
closest point mapping is not parallel to the grid lines. 


n 

Relative Error 

Order 

64 

3.7159 X lO-'" 

- 

128 

2.5786 X 10“^ 

7.17 

256 

4.2361 X 10“'= 

-4.04 

512 

3.2246 X 10“*' 

0.39 

1024 

1.8876 X 10"“ 

0.77 

2048 

1.0132 X 10“^ 

0.90 

4096 

5.2372 X 10-^ 

0.95 

8192 

2.6615 X 10“' 

0.98 


Table 3 Errors for a torus. Relative errors in the numerical approximation of the surface area of a 
torus centered at 0. The distance from the center to the tube that form the torus \s R = 0.75 and the 
radius of the tube is r = 0.25. In this computation, we summed up grid points that are within e = 0.2 
distance from the surface for REco and RE 4 , and e = 0.03 for REi. REoo. RE 4 and REi are the 
relative error using the exact signed distance function, the relative error using a fourth order accurate 
signed distance function and the relative error using a first order accurate signed distance function 
respectively. The Jacobian matrix is approximated by a standard third order accurate differencing 
except for REi where we used a second order accurate differencing to approximate Pp. 


n 

REoo 

Order 

RE 4 

Order 

REi 

Order 

32 

6.2030 X 10“^ 

- 

1.1699 X 10“^ 

- 

5.8000 X 10"^ 

- 

64 

1.8073 X lO”-* 

5.10 

1.0169 X lO-"* 

3.52 

1.4456 X lO""’ 

2.00 

128 

6.6838 X 10-'’ 

4.76 

1.3568 X 10-^ 

6.23 

3.9830 X 10-^ 

1.86 

256 

4.1530 X 10"^ 

4.01 

7.1567 X 10-^ 

4.24 

1.4391 X 10"^ 

1.47 

512 

5.0379 X 10-® 

3.04 

6.1982 X 10““ 

3.53 

5.1463 X lO"'* 

1.48 


Table 4 Errors for a quarter sphere. Relative errors in the numerical approximation of the surface area 
of a quarter sphere with radius R = 0.75 centered at 0. In this computation, we summed up grid 
points that are within e = 0.2 distance from the surface. We used the standard central difference 
scheme to compute each entry of the Jacobian matrix Pp. 


n 

Relative Error 

Order 

32 

9.2825 X lO-'^' 

- 

64 

1.8365 X 10-^ 

2.34 

128 

2.7726 X 10-"^ 

2.73 

256 

7.1886 X 10-^ 

1.95 

512 

1.4811 X lO-'^ 

2.30 


Table 5 Errors for a three quarter sphere. Relative errors in the numerical approximation of the 
surface area of a three quarter sphere with radius R = 0.75 centered at 0 (this is the portion of a 
sphere that misses half of a hemisphere). In this computation, we summed up grid points that are 
within e = 0.2 distance from the surface. Due to this setup, the closest point mapping has a 
discontinuity that stems out from the boundary of the surface. We used the discretization described 
in Section 3.3 to compute each entry of the Jacobian matrix Pp. 


n 

Relative Error 

Order 

32 

1.1726 X 10"^ 

- 

64 

1.1733 X 10-^ 

3.32 

128 

9.1325 X 10-"‘ 

0.36 

256 

3.8238 X 10-"^ 

1.26 

512 

7.8308 X lO-'^ 

2.29 


Table 6 Errors for a coil. Relative errors in the numerical approximation of a coil wrapped around a 
helix. In this computation, we used a constant width for the tubular neighborhood e = 0.1 and took 
the averaging kernels to be defined in (19). 


n 

Relative Error 

Order 

60 

5.5078 X 10-^ 

- 

120 

1.1476 X lO”"* 

2.63 

240 

2.3409 X 10-"^ 

2.29 

480 

3.7166 X lO-'^ 

2.66 





















































